Chapter 3 
Interpolation 



Interpolation is the problem of fitting a smooth curve through a given set of 
points, generally as the graph of a function. It is useful at least in data analy- 
sis (interpolation is a form of regression), industrial design, signal processing 
(digital-to-analog conversion) and in numerical analysis. It is one of those 
important recurring concepts in applied mathematics. In this chapter, we 
will immediately put interpolation to use to formulate high-order quadrature 
and differentiation rules. 

3.1 Polynomial interpolation 

Given + 1 points Xj E M., < j < N , and sample values i/j = f{xj) of 
a function at these points, the polynomial interpolation problem consists in 
finding a polynomial Pn{x) of degree N which reproduces those values: 

In other words the graph of the polynomial should pass through the points 
{xj,yj). A degree-N polynomial can be written as Pn{x) = Xl^Lo^™-^" 
some coefficients ao, . . . ,aN. For interpolation, the number of degrees of 
freedom {N + 1 coefficients) in the polynomial matches the number of points 
where the function should be fit. If the degree of the polynomial is strictly 
less than N, we cannot in general pass it through the points {xj,yj). We 
can still try to pass a polynomial (e.g., a line) in the "best approximate 
manner", but this is a problem in approximation rather than interpolation; 
we will return to it later in the chapter on least-squares. 
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Let us first see how the interpolation problem can be solved numerically 
in a direct way. Use the expression oi pn into the interpolating equations 

yj=PNixjy. 

N 

^anx]^yj, j^O,...,N. 

n=0 

In these + 1 equations indexed by j, the unknowns are the coefficients 
ao, . . . , Cat. We are in presence of a linear system 

N 

Va = y, ^ ^ Vjnan = 

n=0 

with V the so-called Vandermonde matrix, Vjn = x^, i.e., 

/l Xq ■ ■ ■ Xq\ 

1 Xi ■ ■ ■ x\ 
. . . . 

\l Xn ■■■ x%J 

We can then use a numerical software like Matlab to construct the vector of 
abscissas xj, the right-hand-side of values yj, the V matrix, and numerically 
solve the system with an instruction like a = \^\y (in Matlab). This gives 
us the coefficients of the desired polynomial. The polynomial can now be 
plotted in between the grid points Xj (on a finer grid), in order to display 
the interpolant. 

Historically, mathematicians such as Lagrange and Newton did not have 

access to computers to display interpolants, so they found explicit (and el- 
egant) formulas for the coefficients of the interpolation polynomial. It not 
only simplified computations for them, but also allowed them to understand 
the error of polynomial interpolation, i.e., the difference f{x) — pn{x). Let 
us spend a bit of time retracing their steps. (They were concerned with 
applications such as ffiting curves to celestial trajectories.) 

We'll define the interpolation error from the uniform {L°°) norm of the 
difference f — Pn'- 

11/ -PjvIIoo := max|/(x) -pjv(x)|, 

X 

where the maximum is taken over the interval [xq^xn]- 
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Call Pn the space of real- valued degree- polynomials: 

AT 



n=0 

Lagrange's solution to the problem of polynomial interpolation is based 
on the following construction. 

Lemma 1. (Lagrange elementary polynomials) Let {xj,j — 0, . . . , N} be a 

collection of disjoint numbers. For each k — 0, . . . , N , there exists a unique 
degree-N polynomial Lk{x) such that 

l^k[Xj) = djk = j Q ifj^k. 
Proof. Fix k. L^ has roots at Xj for j ^ k, so Lk must be of the form^ 



Lk{x) = C JJ(a; - Xj). 
Evaluating this expression at x — Xk, we get 

Hence the only possible expression for Lk is 



Lk{x) = 



These elementary polynomials form a basis (in the sense of linear algebra) 
for expanding any polynomial interpolant p^. 

^That's because, if we fix j, we can divide Lk{x) by {x — Xj), j ^ k. We obtain 
Lk[x) = (x — Xj)q(x) + r(x), where r(x) is a remainder of lower order than x — Xj, i.e., 
a constant. Since Li^{xj) = we must have r{x) = 0. Hence {x — Xj) must be a factor 
of Lk{x). The same is true of any {x — Xj) for j ^ k. There are N such factors, which 
exhausts the degree A''. The only remaining degree of freedom in is the multipUcative 
constant. 
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Theorem 4. (Lagrange interpolation theorem) 

Let {xj,j = 0, . . . , A^} be a collection of disjoint real numbers. Let {yj,j = 
0, . . . , A^} be a collection of real numbers. Then there exists a unique G -Pjv 
such that 

Its expression is 

N 

Pn{x) = ^ykLk{x), (3.1) 

k=0 

where Lk{x) are the Lagrange elementary polynomials. 
Proof. The justification that (3.1) interpolates is obvious: 

N N 

PN{xj) = ^ykLk{xj) ^^ykLk{xj) = yj. 

k=0 k=0 

It remains to see that Pat is the unique interpolating polynomial. For this 
purpose, assume that both pjy and gjv take on the value yj at Xj. Then 
tn = Pn — In is a polynomial of degree N that has a root at each of the 
A'" + 1 points xo, . . . , xjv- The fundamental theorem of algebra, however, says 
that a nonzero polynomial of degree N can only have N (complex) roots. 
Therefore, the only way for r^r to have A'^ + 1 roots is that it is the zero 
polynomial. So p^ = q^. 



By definition, 

N 

Pn{x) ^^f{xk)Lk{x) 

k=0 

is called the Lagrange interpolation polynomial of / at xj. 
Example 7. Linear interpolation through {xi,yi) and {x2,y2)-' 

Li{x) = , L2{x) 



Xi -X2 X2- Xi 



Pi{x) = y\Li{x) + y2L2{x) 



y2-yi , yiX2 - y2Xi 

X H 

X2 - Xi X2 - Xi 

, 1/2 - yi . V 

2/1 H [x - xi). 

X2 - Xi 
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Example 8. (Example (6.1) in Suli-Mayers) Consider f{x) = , and inter- 
polate it by a parabola {N — 2) from three samples at xq = —1, xi — 0,X2 — 1. 
We build 

(x-xi){x-X2) 1 , . 

Lq(x] — 7^ rj — -Xix — 1) 

^ ^ {Xo - Xi){Xo - X2) 2 ^ ^ 

Similarly, 

Li(x) = l-x^, L2{x) ^ -x{x + 1). 
So the quadratic interpolant is 

P2{x) = e-^Loix) + e°Li(a;) + e^L2(x), 
= 1 + sinh(l) X + (cosh(l) — 1) x^, 
~ 1 + 1.1752 x + 0.5431 

Another polynomial that approximates e^ reasonably well on [—1, 1] is the 
Taylor expansion about x — 0: 

t2{x) = 1 + X+ y. 

Manifestly, p2 is not very different from t2. 
(Insert picture here) 

Let us now move on to the main result concerning the interpolation error 
of smooth functions. 

Theorem 5. Let f e C^+^[a, b] for some N > 0, and let {xj : j = 0,. . . ,N} 
be a collection of disjoint reals in [a, b]. Consider p^ the Lagrange interpola- 
tion polynomial of f at Xj. Then for every x e [a, b] there exists ^{x) e [a, b] 
such that 

f{x) -Pn{x) = (AT + l)! 

where 

N+l 

nN+l{x) =W{X-Xj). 

i=i 

An estimate on the interpolation error follows directly from this theorem. 

Set 

M^+i= max 1/(^+1) (x) I 
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(which is well defined since is continuous by assumption, hence reaches 

its lower and upper bounds.) Then 

\fix) -Pn{x)\ < J^^^\^N+i{x)\ 

In particular, we see that the interpolation error is zero when x — Xj for 
some j, as it should be. 

Let us now prove the theorem. 

Proof. (Can be found in Suli-Mayers, Chapter 6) 

In conclusion, the interpolation error: 

• depends on the smoothness of / via the high-order derivative /^^"^^^; 

• has a factor l/(iV +1)! that decays fast as the order iV — > cxd; 

• and is directly proportional to the value of ttn^iIx), indicating that the 
intcrpolant may be better in some places than others. 

The natural follow-up question is that of convergence: can we always 
expect convergence of the polynomial interpolant as N ^ oo? In other 
words, does the factor 1/(A'" -|- 1)! always win over the other two factors? 

Unfortunately, the answer is no in general. There are examples of very 
smooth (analytic) functions for which polynomial interpolation diverges, par- 
ticularly so near the boundaries of the interplation interval. This behavior 
is called the Runge phenomenon, and is usually illustrated by means of the 
following example. 

Example 9. (Runge phenomenon) Let f{x) for x G [—5,5]. Interpolate it 
at equispaced points Xj = lOj/N, where j = —N/2, . . . , N/2 and N is even. 
It is easy to check numerically that the interpolant diverges near the edges of 
[—5, 5], as N —>■ oo. See the Trefethen teoctbook on page 44 for an illustration 
of the Runge phenomenon. 

(Figure here) 

If we had done the same numerical experiment for a; G [— 1, 1], the inter- 
polant would have converged. This shows that the size of the interval matters. 
Intuitively, there is divergence when the size of the interval is larger than the 
"features" , or characteristic length scale, of the function (here the width of 
the bump near the origin.) 
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The analytical reason for the divergence in the example above is due in 
no small part to the very large values taken on by 7[n^i{x) far away from the 
origin — in contrast to the relatively small values it takes on near the origin. 
This is a problem intrinsic to equispaced grids. We will be more quantitative 
about this issue in the section on Chebyshev interpolants, where a remedy 
involving non-equispaced grid points will be explained. 

As a conclusion, polynomial interpolants can be good for small N, and 
on small intervals, but may fail to converge (quite dramatically) when the 
interpolation interval is large. 



3.2 Polynomial rules for integration 

In this section, we return to the problem of approximating J^u(x)dx by 
a weighted sum of samples u{xj), also called a quadrature. The plan is 
to form interpolants of the data, integrate those interpolants, and deduce 
corresponding quadrature formulas. We can formulate rules of arbitrarily 
high order this way, although in practice we almost never go beyond order 4 
with polynomial rules. 



3.2.1 Polynomial rules 

Without loss of generahty, consider the local interpolants of u{x) formed near 
the origin, with xq — 0,xi — h and X-i — —h. The rectangle rule does not 
belong in this section: it is not formed from an interpolant. 

• The trapezoidal rule, where we approximate u{x) by a line joining 
(0, w(xo)) and {h,u{xi)) in [0, /i]. We need 2 derivatives to control the 
error: 

u[x) — pi[x) H x[x — n), 

= u{^)Lq{x) + u{h)Li{x), 

ph ph 

/ Lo(x)dx = / Li[x)dx = h/2, (areas of triangles) 
Jo Jo 

r I'^^^^^xix -h)\dx< Cmax|M"(C)|/i'. 
Jo ^ ^ 
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The result is 



As we have seen, the terms combine as 



1 h ^"^ h 

u{x) dx — —u{xq) + h u{xj-\) + —u{xn) + 0{h^). 



Simpson's rule, where we approximate u{x) by a parabola through 
{—h,u{x_i)), {0,u{xo)), and {h,u{xi)) in [—h,h]. We need three deriva- 
tives to control the error: 

u(x) = p2(x) + ^ ^i^^''\ x + h)x(x - h), 
b 

P2{x) = u{-h)L_i{x) + u{0)Lo{x) + u{h)Li{x), 

r , , xix-h) T / \ (x + h)(x-h) T / \ (x + h)x 
L-i(x) = Lo(x) = — , L,(x) = 

/h rh rh 

L^i(x)dx — / Li(x)dx — h/3, / Lo(x)dx — Ah/S, 
-h J-h J-h 

'\u"'{i{x)) 



Jo 



'0 

The result is 



g ,x + h)x(x -h)\dx< Cmax |ii'"(OI^ ■ 



The composite Simpson's rule is obtained by using this approximation 
on [0, 2h] U [2h, Ah]Li . . .[1 ~ 2h, 1], adding the terms, and recognizing 
that the samples at 2nh (except and 1) are represented twice. 

^ / X , , fuiO) + 4:u(h) + 2u(2h) + Au(3h) + ... + 2u(l - 2h) + Auil - h) + u(l) 
u[x)dx — h ' 

It turns out that the error is in fact 0{h^) on [— /i, /i], and 0{h^) on 
[0, 1], a result that can be derived by using symmetry considerations 
(or canceling the terms from Taylor expansions in a tedious way.) For 
this, we need u to be four times differentiable (the constant in front of 
involves u"" .) 
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The higher-order variants of polynomial rules, called Ncwton-Cotes rules, 
are not very interesting because the Runge phenomenon kicks in again. Also, 
the weights (like 1,4,2,4,2, etc.) become negative, which leads to unacceptable 
error magnification if the samples of u are not known exactly. 

Piecewise spline interpolation is not a good choice for numerical integra- 
tion either, because of the two leftover degrees of freedom, whose arbitrary 
choice affects the accuracy of quadrature in an unacceptable manner. (We'll 
study splines in a later section.) 

We'll return to (useful!) integration rules of arbitrarily high order in the 
scope of spectral methods. 



3.3 Polynomial rules for differentiation 

A systematic way of deriving finite difference formulas of higher order is to 
view them as derivatives of a polynomial intcrpolant passing through a small 
number of points neighboring xj. For instance (again we use —h,0,h as 
reference points without loss of generality): 

• The forward difference at is obtained from the line joining (0,m(0)) 
and {h, u{h))\ 

pi{x) = u{Q))Lq{x) +u{h)Li{x), 
Lo{x) = -J^^ -^i(^) = 

We already know that u'{0) — p'i(O) — 0{h). 

• The centered difference at is obtained from the line joining {—h, u{—h)) 
and {h,u{h)) (a simple exercise), but it is also obtained from differen- 
tiating the parabola passing through the points {—h,u{—h)), (0,ti(0)), 
and {h,u{h)). Indeed, 

P2{x) = u{—h)L_i{x) + u{0)Lo{x) + u{h)L-i{x), 

^ , , x(x — h) ^ , , (x + h)(x — h) ^ , , (x + h)x 
L-A^) = L,{x) = , L,{x) = 
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We already know that u'{0) -^2(0) = 0{h'^). 

• Other examples can be considered, such as the centered second differ- 
ence (-1 2 -1), the one-sided first difference (-3 4 -1), etc. 

Differentiating one-sided interpolation polynomials is a good way to ob- 
tain one-sided difference formulas, which comes in handy at the boundaries 
of the computational domain. The following result establishes that the order 
of the finite difference formula matches the order of the polynomial being 
differentiated. 

Theorem 6. Let f G C'^'^^[a,b], {xj,j = 0,...,N} some disjoint points, 
and pjv the corresponding interpolation polynomial. Then 

for some ^ G [a, b], and where 7r^(x) = {x — rji) . . . {x — tjn) for some rjj e 



\x 



^l,Xj]. 



The proof of this result is an application of Rolle's theorem that we leave 
out. (It is not simply a matter of differentiating the error formula for poly- 
nomial interpolation, because we have no guarantee on d^/dx.) 

A consequence of this theorem is that the error in computing the deriva- 
tive is a 0{h^) (which comes from a bound on the product 7r^(a;).) 

It is interesting to notice, at least empirically, that the Runge's phe- 
nomenon is absent when the derivative is evaluated at the center of the 
interval over which the interpolant is built. 



3.4 Piecewise polynomial interpolation 

The idea of picccwisc polynomial interpolation, also called spline interpola- 
tion, is to subdivide the interval [a, b] into a large number of subintervals 
[xj_i,Xj], and to use low-degree polynomials over each subintervals. This 
helps avoiding the Runge phenomenon. The price to pay is that the inter- 
polant is no longer a function - instead, we lose differentiability at the 
junctions between subintervals, where the polynomials are made to match. 
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If we use polynomials of order n, then the interpolant is picccwise C°°, 
and overall C^, with k < n — 1. We could not expect to have k = n, because 
it would imply that the polynomials are identical over each subinterval. We'll 
see two examples: linear splines when n — 1, and cubic splines when n — 3. 
(We have seen in the homework why it is a bad idea to choose n — 2.) 

3.4.1 Linear splines 

We wish to interpolate a continuous function f{x) of a; G [a,b], from the 
knowledge of f{xj) at some points Xj, j — 0, . . . ,N, not necessarily equis- 
paced. Assume that xq = a and xn = b. The piecewise linear interpolant 
is build by tracing a straight line between the points (xj-i, f{xj-i)) and 
{xj, f{xj))); for j = 1, . . . , the formula is simply 

Sl{x) = ^^'^ ^ f{Xj-l) + ^ f{Xj), X e [Xj-i, Xj]. 

X j X j ]^ X j X j ]_ 

In this case we see that sl{x) is a continuous function of x, but that its 
derivative is not continuous at the junction points, or nodes xj. 

If the function is at least twice differentiable, then piecewise linear inter- 
polation has second-order accuracy, i.e., an error 

Theorem 7. Let f e C^[a, b], and let h = maXj=i^,,_^N{xj — Xj^i) be the grid 
diameter. Then 

h^ 

\\f - SL\\L--[a,b] < yll/'Ik-M- 

Proof. Let x e [xj_i,Xj] for some j = 1, . . . , N. We can apply the basic result 

of accuracy of polynomial interpolation with n — 1: there exists ^ e [ X j ]_ ^ X j J 

such that 

f{x) - Sl{x) = ^/"(O (X - Xj-i){x - Xj), X e [Xj-i,Xj]. 

Let hj = Xj — Xj^i. It is easy to check that the product {x — — Xj) 

takes its maximum value at the midpoint ^^']^^^ ^ and that the value is h'j/A. 
We then have 

\f{x) - SLix)\ < -f max |/"(0I, x e [xj-i,Xj]. 

The conclusion follows by taking a maximum over j. 
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We can express any piecewise linear interpolant as a superposition of 
"tent" basis functions 4>k{x)- 

k 

where (pki^) is the piecewise hnear function equal to 1 at a; = x^, and equal 
to zero at all the other grid points Xj, j ^ k. Or in short, (f)k{xj) — 5jk- On 
an equispaced grid xj — jh, an explicit formula is 

where 

S(i){x) = x+- 2{x - h)+ + {x - 2h)+, 

and where x+ denotes the positive part of x (i.e., a; if a; > 0, and zero if 
X < 0.) 

Observe that we can simply take = f{xk) above. The situation will be 
more complicated when we pass to higher-order polynomials. 

3.4.2 Cubic splines 

Let us now consider the case n = 3 of an interpolant which is a third-order 
polynomial, i.e., a cubic, on each subinterval. The most we can ask is that the 
value of the interpolant, its derivative, and its second derivative be continuous 
at the junction points xj. 

Definition 5. (Interpolating cubic spline) Let f G C[a,b], and {xj,j = 
0, . . . , N} C [a, b]. An interpolating cubic spline is a function s{x) such that 

L s{xj) = f{xj); 

2. s{x) is a polynomial of degree 3 over each segment [xj_i,Xj\; 

3. s{x) is globally C^, i.e., at each junction point Xj, we have the relations 

s{x-) = s{xl), s'ixj) = s'ixl), s"{x-) = s"{xl), 

where the notations s{xj) and s{Xj') refer to the adequate limits on the 
left and on the right. 
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Let us count the degrees of freedom. A cubic polynomial has 4 coefficients. 
There are + 1 points, hence N subintervals, for a total of AN numbers to 
be specified. 

The interpolating conditions s{xj) = f{xj) specify two degrees of freedom 
per polynomial: one value at the left endpoint Xj^i, and one value at the right 
endpoint xj. That's 2N conditions. Continuity oi s{x) follows automatically. 
The continuity conditions on s', s" are imposed only at the interior grid points 
Xj for j = 1, . . . , — 1, so this gives rise to 2(A^ — 1) additional conditions, 
for a total of AN — 2 equations. 

There is a mismatch between the number of unknowns (4A") and the 
number of conditions (4A^ — 2). Two more degrees of freedom arc required to 
completely specify a cubic spline interpolant, hence the precaution to write 
"an" interpolant in the definition above, and not "the" interpolant. 

The most widespread choices for fixing these two degrees of freedom are: 

• Natural splines: 

s"{xo) = s"{xn) = 0. 

If s{x) measures the displacement of a beam, a condition of vanishing 
second derivative corresponds to a free end. 

• Clamped spline: 

s'{xo)=po, s'{xn)^Pn, 

where Pq and pjv are specified values that depend on the particular 
application. If s{x) measures the displacement of a beam, a condition 
of vanishing derivative corresponds to a horizontal clamped end. 

• Periodic spline: assuming that s{xo) — s{xn), then we also impose 

s'{xo) = s'{xn), s"{xo) = s"{xn)- 

Let us now explain the algorithm most often used for determining a cubic 
spline interpolant, i.e., the 4A^ coefficients of the A^ cubic polynomials, from 
the knowledge of f{xj). Let us consider the natural sphne. 

It is advantageous to write a system of equations for the second derivatives 
at the grid points, that we denote aj = s"{xj). Because s{x) is piecewise 
cubic, we know that s"{x) is piecewise linear. Let hj = Xj — Xj^i. We can 
write 

s {x) = aj-i + aj, X e [xj-uXj]. 
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Hence 



S[X) 



Gj-i + — + aj{x - Xj-i) + (5j{xj - x). 



6hj ' 6hj 



We could have written ax + b for the effect of the integration constants in 
the equation above, but writing it in terms of aj and (3j makes the algebra 
that follows simpler. The two interpolation conditions for [ be 
written 



2 



aj-ih 

s{xj-i) = f{xj-i) f{xj-i) = — - — + hj(3j, 

= f{xj) ^ f{xj) = + hjaj. 

One then isolates aj, j3j in this equation; substitutes those values in the 
equation for s{x)] and evaluates the relation s'{x^) = s'{xj'). Given that 
Co = cTAf = 0, we end up with a system of — 1 equations in the N — 1 
unknowns ai, . . . ,aN- Skipping the algebra, the end result is 

hjaj^i + 2{hj+i+hj)aj + hj+iaj+i = 6 ^ j— — j 

We are in presence of a tridiagonal system for aj. It can be solved efficiently 
with Gaussian elimination, yielding a LU decomposition with bidiagonal 
factors. Unlike in the case of linear splines, there is no way around the fact 
that a linear system needs to be solved. 

Notice that the tridiagonal matrix of the system above is diagonally dom- 
inant (each diagonal element is strictly greater than the sum of the other 
elements on the same row, in absolute value), hence it is always invertible. 

One can check that the interpolation for cubic splines is 0{h'^) well away 
from the endpoints. This requires an analysis that is too involved for the 
present set of notes. 

Finally, like in the linear case, let us consider the question of expanding 
a cubic spline interpolant as a superposition of basis functions 

k 

There are many ways of choosing the functions (f)k{x), so let us specify that 
they should have as small a support as possible, that they should have the 
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same smoothness as s{x) itself, and that they should be translates of each 
other when the grid xj is equispaced with spacing h. 

The only solution to this problem is the cubic B-spline. Around any 
interior grid point x^, is supported in the interval [xk-2,Xk+2], and is given 
by the formula 

where 

5(3)(x) ^xl- A(x - h)l + 6(x - 2h)l - A(x - 3h)l + (x - Ah)l, 

and where a;+ is the positive part of x. One can check that (j)k{x) takes the 
value 1 at Xk, 1/4 at Xk±i, zero outside of [xk-2,Xk+2], and is at each 
junction. It is a bell-shaped curve, it is an interesting exercise to check that 
it can be obtained as the convolution of two "tent" basis functions of linear 
interpolation: 

= cS(i){x) * S(^i){x), 

where c is some constant, and the symbol * denotes convolution: 

f*g{x) = / f{y)g{x-y)dy. 
Jr 

Now with cubic B-splines, we cannot put Cfc = fi^k) anymore, since 
(f)k{xj) 7^ Sjk- The particular values of Ck arc the result of solving a lin- 
ear system, as mentioned above. Again, there is no way around solving a 
linear system. In particular, if f{xk) changes at one point, or if we add a 
datum point {xk, fk), then the update requires re-computing the whole inter- 
polant. Changing one point has ripple effects throughout the whole interval 
[a, b]. This behavior is ultimately due to the more significant overlap between 
neighboring (jik in the cubic case than in the linear case. 
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